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Abstract 

Structured materials like photonic crystals require for optimal use a high precision both 
on position and optical characteristics of the components which they are made of. Here, we 
present a simple tomographic algorithm, based on a specific Green's function together with a 
first-order Born approximation, which enables us to localize and characterize identical defects 
in finite-size photonic crystals. This algorithm is proposed as a first step to the monitoring of 
such materials. Illustrative numerical results show in particular some possibility of focalization 
beyond the Rayleigh criterion. 

1 Introduction 

Photonic crystals (PCs) [HE] are periodic, dielectric or metallic structures, which possess a variety 
of band dispersions and band gaps. They are found in antennas [31 H], waveguides [5], negative 
refractive index materials [6] to quote a few usages. As is well-known, such structures require for 
optimal use that both position and optical properties of materials they are consisting of be very 
precise. Though the density of states is zero within the photonic band gap, by perturbing a single 
lattice site, a single mode or a set of closely spaced modes which have frequencies within this gap 
is permitted, e.g., a single column can be removed from the crystals or replaced with another the 
size, shape, or dielectric constant of which is different from the original. Properties [HI El E] an d 
modeling PUIE] of PCs have been extensively studied, though mostly for 2-D configurations, these 
last two decades. 

The present investigation is intended to be a first step towards the health monitoring of such 
structures. A low-complexity tomographic algorithm that enables us to localize and characterize 
simple defects consisting in either an absence of circular cylinders or a modification of the optical 
properties of such cylinders in finite-sized 2-D photonic crystals. This is done from a low-frequency 
excitation outside the usual frequency range of the band gap. Both the specific Green's function and 
the response to the low-frequency excitation of the original, intact structure (and of the damaged 
structure) can be calculated by means of the so-called multipole method [T2], and they can be 
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used to solve the inverse problem via a first-order Born approximation. The introduction of a 
specific Green's function directly exhibits the discrepancies between the initial configuration and 
the configuration with defects |13j . 

In the background of this investigation, recent studies have also shown that the utilization of 
a structured embedding medium, like PCs, could lead to focusing and resolution of a tomographic 
inversion algorithm beyond the Rayleigh criterion [14|, I15j . The possibility to encounter such a 
phenomenon is investigated in addition herein. 



2 Field (Data) and specific Green's function. 

Let us consider a structure made of ./V parallel circular cylinders C 3 , identified by superscript j £ N, 
of radius R? and of optical index rf , located at r 3 = (r 3 , 9 3 ) in the global polar coordinate system 
in the cross-sectional plane (this is a 2-D scattering configuration). 

As indicated in the above, both fields and a specific Green's function (the latter being the field 
solution of the problem when the structure is excited by a given line source) are calculated by means 
of the multipole method [12]. Key to this approach are the local field expansions or multipole 
expansions in the vicinity of each cylinder in the polar coordinate system linked to that cylinder 
(which are derived from the application of the Graf's addition theorem |16j): 



E e {v x )= [B l m H£Hkn) + A l m J m (kn)]e ime < (1) 

m=— oo 

wherein H$ is the first-kind Hankel function of order m, J m is the Bessel function of order m, 
B l m are the coefficients of the scattered field by the Z-th cylinder, A l m are those of the incident 
field impinging upon the Z-th cylinder, and ri = (rj,fy) are the coordinates of a point P in the 
polar coordinate system linked to this l-th. cylinder. The local incident field on the l-th cylinder is 
generated by the actual incident field E inc as well as by the fields that are scattered by all other 
cylinders j, j ^ I. Their coefficients also take the form 

TV oo 

= Km + ^mpBp (^) 

j=l,j^lp=-oo 

wherein K l m are the coefficients of the actual incident field (either a planar incident wave in the cross- 
sectional plane, K l m = (— i) m exp (— ifcr'cos (# mc — — im8 mc ^, or a cylindrical wave generated 

by an exterior line source set parallel to the axis of the cylinder) and Sm P = H^_ p (^krj^j e^ p ^ m ^ 8 ' 

are translation terms, (rj,9j) being the coordinates of the j-th cylinder in the polar coordinate 
system associated to the l-th cylinder. 

Coefficients A l m and B l m are related together via the continuity of the tangential components 
of the electric and magnetic fields to be imposed at the cylinder boundaries. To derive these 
relationships, the interior field expansion within the cylinder I is used as 

oo 

= [0LHg> (krfn) + C l m J m (k V l n)] e im6 ' (3) 
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wherein C l m are the coefficients of the scattered field inside the I cylinder and Q l m = jx l Jm(kr] l rf)e ime i 
are the coefficients of a field generated by a line source located at (rf,0f) inside the cylinder / in 
the polar coordinate system associated to it. The presence of the interior source is indicated by 
the term \ l valued to 1 when the source is present and otherwise. 

The continuity conditions at the boundaries are most conveniently expressed in terms of cylin- 
drical harmonic reflection and transmission coefficients [7] as 

J}1 T>1 Al i rpl r\\ 

m - n 'm m ' m^m ( a\ 

C l = T' 1 A 1 4- T}' 1 f) 1 ^ ' 

m m ' ^m^m' 

In the above, vectors B = [B l m ], K = [K l m ] and Q = [Q l m ], as well as matrices S = [S^tp], 
R = diagit!^ (with identical definition applying to the other reflection and transmission matrices 
and the other coefficients), are introduced in order to deduce from ([2]) and (jl]) the system of linear 
equations in the source coefficients B as 

(I - RS) B = RK + TQ. (5) 

Upon solving the above linear system both in the interior (using the second equation of @) 
and in the exterior, the specific Green's function is made available. 

3 A simple tomographic algorithm. 

From now on, one is considering a finite-size crystal (FSC) with hexagonal symmetry as is sketched 
in Fig. [H It is made of iV (iV will be chosen as 85 in the numerical examples) circular cylinders of 
same radius R and same optical index rj. The cylinders C-?, j G [lj-^V] are ordered such that C 1 is 
located at the bottom left corner and C N at the top right corner. The distance between the centers 
of the closest cylinders is denoted by d. The electric field is calculated on a circle of radius r = 20d 
(this is a rather arbitrary value, what matters is that one stays fully outside the crystal). Fields 
are time-harmonic, with wavelength in free space as A, and wavenumber k = 2ir/X. A defect is 
obtained either by removing the Z-th cylinder, or by modifying its optical index, de facto creating 
a novel configuration (denoted as Z-th DFSC). 




Figure 1: Finite crystals of N = 85 parallel circular cylinders with a single defect. The circle 
around the structure is the one used for the field computation. 
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Tomographic algorithms are usually derived from an integral formulation of the solution based 
on the so-called background (free- field) Green's function. Use of a specific Green's function, i.e., 
the solution (for a given line source) of a problem close to the one at hand reduces the kernel of 
the integral. Here, since the structure of materials which one is interested in is designed to exhibit 
specific properties, it is consistent to consider the Green's function of the configuration associated 
with this structure (i.e., in the absence of defects) whenever intending to carry out the health 
monitoring. 



3.1 Localization and characterization of a single defect 

One introduces g(r,rf) and E(r) as the specific Green's function and the total electric field 
calculated for the FSC, and E l (r) the total electric field calculated for the l-th DFSC. These 
fields are related by the well-known Fredholm integral equation of first kind E l (r) — E(r) = 

J cl g(r,Tf)k 2 ([ V l ) 2 -r,A E l (vf)dS. 



Then, a first-order Born approximation, whose accuracy is related to the fact that the field at the 
location of the defect should be weakly modified by the defect, is employed. Any mode associated 
to the defect, possibly leading to entrapment of the field, should be excited. Complementarily, 
one applies a low-frequency approximation, which requires that the radius of each cylinder is small 
enough with respect to the source wavelength. Furthermore, one works outside the band gap. 
Under these conditions the Green's function and the field E l (rf) read as 

lJ r„ tU~*\X^lj (J\ rr(l) 



g(r,r\) = -fj (fatf) £ 6& (x l ) (fc|r - H| y ^ 

E l (vl) = 4 Jo (fe/rf) =TvA l Jo (kwl) 

In Eq.© one has introduced 6q, such that \TqJq {krjrf) 6q = Bq is satisfying the linear set of 
equations (I — RS) b = with vector x = [x 1 ]- Coefficients of the locally impinging field A l can be 
calculated via the solution of Eq.@ with B{ satisfying (I - RS) B = R7C, K& = e~' lkrj «*(« <no -« i ). 

By making use of relation [17J J x (Jo(ax)) 2 dc = ^- |(Jo(ax)) 2 + (Ji(ax)) 2 |, the first-order 
Born approximation of the integral equation is 

N 

E\v) - E(r) » V( l A l £ (x l ) (k\v - r»|) (7) 

where V = iT T' Q (kR) 2 {( J {k7]R)) 2 + ( Ji(fcr/i?)) 2 } /4 and C' = (Vf - ?? 2 )- V is found to be 
independent from both the location and the material characteristics of the defect, while £ is the 
contrast function. Localization and characterization of the defect can thus be fully decoupled. 

Let us introduce the normalized vector v = V/V • V* (V* being the complex conjugate of 
V) associated with vector V = \(E l (r) — E(r)} and iV normalized vectors gj = Gj/Gj ■ Gj 



associated with vector G J 



The defect is localized whenever 



= 1 1/(1 — ||^ J '||)|, letting zi = gj • v*/gj ■ gj , is maximum. The function is derived from 
The parameter corresponds with focalization at the defect location. 
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Finally, the optical index of the l-th cylinder can be retrieved (this works quite well as seen 
next, yet it remains heuristic) by averaging the value of (V(r)/G l {r) + rj 2 ) 2 over the measured 
data, i.e. rf = mean ( (V(r)/G l (r) + ?y 2 ) 2 ), wherein mean(x(r)) means the average value of x(r) 



over the measured data. 

3.2 Localization and characterisation of two identical defects 

The same assumptions and procedure as above is followed for two identical defects, which corre- 
sponds to the so-denoted (i,j)-th DFSC configuration. The first-order Born approximation now 
requires that the field at the location of one defect be only weakly modified by this defect but also 
by the other one. This means that the defects are "well separated", i.e., they are not interacting 
together. The approximation of the integral equation becomes 

N 

E^\r) - E(r) « 2>c£ (x*) + 4tf> U')) ^ (*k " r 'l) ( 8 ) 

1=1 

Similarly with what has been done in the previous subsection, one is introducing the normalized 
vector v associated with vector V = UE^'ir) — E(r)) /P] and N x N normalized vectors g( q '') 

associated with vector G^' 1 ) = \j2p=i {A b l (x q ) + -4^0 [x 1 )) (k\r - r p |)j . The defect is now 

localized whenever V^ q ^ = |1/(1 — ||z^'^||)|, letting z^ q ^ = g( q '') • v*/g( q ' 1 ) • g^ 01 ' 1 )*, is maximum. 

The optical index of the (z,j)-th cylinders then follows as rf^ = mean I (V (r) / G^'^ (r) + r/ 2 ) 3 J. 



4 Numerical results 

Data are computed by use of the multipole method. The infinite sums X^m=-oo are truncated to 
^2m=-M sucn that M = int (^(krjR)~3 + kr]R + 5^ , where int(x) is the entire part of x. For the 
inverse problem, one does not make use of the fact that the scattered field is isotropic for the calcu- 
lation of the specific Green's function as it is proposed in [12] in the low-frequency approximation. 
The latter is calculated directly, by multiplying the vector b by a matrix H = [H^ (k\r — r J |)] ; 
which is stored once. In the same fashion, the available asymptotic formulae of the reflection and 
transmission coefficients are not employed; the formulae of these coefficients are as in [7]. 

The wavelength and the optical index ry are set to A = 20 and r\ = 2.9. The low- frequency 
approximation is then valid for R small enough to have kr/R « 1. In the following one assumes 
that R = 0.15 and one mostly investigates the localization and characterization of defects located 
in the central part of the PCs when tt/2. 
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(a) (b) 




Figure 2: Crystal with one defect: R = 0.15, r/ 43 = 1 and 6 inc = vr/2. (a) V j x 10~ 3 when d = 4 
(^ 43 = 1.1) and (b) V j x 10~ 3 when d = 1 (^ 43 = 1.08). 

Figure [2] shows "P 5 for a 43-th DFSC, with rf = 1 (i.e., the 43-th cylinder is removed) and 
inc = vr/2 when d = 4 (i?/d = 3.75 x 10~ 2 ) and d = 1 {R/d = 0.15). In both cases, the modified 
cylinder is clearly retrieved and rf is found with a relative error on its real part Er = (jR(rj l ) — rf) /rf 
less than 0.1. The 43-th cylinder seems (in the sense that the function V 1 does not point to another 
cylinder) to be retrieved with an accuracy of d that is much smaller than the Rayleigh criterion 
A/2 = 10. 

However, this is partly due to the representation chosen, i.e., the value of V is very large 
since ||^|| ~ 1, focusing beyond the Rayleigh criterion playing its part only to some extent. Indeed, 
strictly speaking, super-resolution would mean that the width at half- height of \\z J \\ is (significantly) 
smaller than half a wavelength. A cut of \\z^\\ along the axis passing through r 43 with an angle 
of vr/6 is displayed in Fig. [3] for 6 inc = vr/2, when d = 2 (R/d = 7.5 x 10~ 2 ) and d = 1. Super- 
resolution, in terms of focusing accuracy, is validated in case d = 1 and not so much in case d = 2. 




|rJ-r a3 | 

Figure 3: Crystal with a single defect: R = 0.15, r/ 43 = 1 and 8 mc = tt/2. Cut of along the 
axis going through r 43 with an angle 9 = tt/6. 

Let us notice that the problem at hand can be interpreted as the retrieval of the location of an 
induced line source within a cylinder, whilst the problem attacked in [TO] consists in the retrieval 
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of a line source located outside all cylinders. In this sense, the problem here is quite different but 
like effects are observed. 




Figure 4: Crystal with a single defect when d = 4: R = 0.15. (a) x 10 3 when r/ 43 = 1 and 

Qinc _ ^4 (~43 _ L1 ) and ( b ) -pj x 1Q -3 when ^77 = j and Qinc = ^ ^7 = 

The retrieved value of rf when d = 4 is the same when 9 mc is varied within [0;7r/2] as well as 
for a defect which is not located in the central part of the PCs, whilst the height of the peak of V 1 
depends on both 9 mc and defect location, Fig. HI with no obvious rule however. 




1 2 3 4 5 

rf 



Figure 5: Crystal with a single defect when d = 4: R = 0.15 and 9 inc = ir/2. Relative error on the 
reconstructed value rf 13 for various r/ 43 . The modified cylinder is systematically retrieved. 

Figure shows the relative error on the reconstruction of r/ 43 for various values of ry 43 , with 
Qinc _ 75-/2, when d = 4. The smaller the contrast is, the better the retrieval of rf is. In 
particular, one is able to retrieve a ry 43 = 2.8, that represents a variation of 3.45% from rj, with a 
relative error less than 10 . The results remain accurate when the low- frequency approximation 
is not anymore valid at the defect location yet remains valid for cylinders that are constituting the 
structured background. 
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Figure 6: Crystal with a single defect when d = 4: rf = 1 and 6 mc = tt/2. Evolution of Er (dashed 
line) and of the height of the peak of V 1 normalized by its value for R = 0.075 (solid line) for 
various values of R. 



The relative height h l of the peak of V 1 is defined by h! = V 1 — mean^ \P^) > wherein mean.,- 



is the average value of x- 7 over j € [1,-/V]. The evolution of Er and of h l with R £ [7.5 x 10 _2 ,0.6] 
is shown in Fig. [6] with 8 mc = tt/2 when d = 4 and r] = 1. It is observed that, even though the 
low-frequency approximation does not hold for the background, the modified cylinder is retrieved, 
but the reconstructed values rf are not accurate anymore. When rj is varied, the reconstructed 
values (in particular the imaginary part) of rf are now inaccurate even for kr]R = 0.27. This means 
that the result, and as a matter of fact the use of the first-order Born approximation, remains much 
more appropriate when R becomes large than when rj becomes large. 
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Figure 7: Crystal with a single defect when d = 4: R = 0.15, rf = 1 and 9 mc = tt/2. (a) the solid 
curve depicts ||v|| and the dashed curve depicts ||v no j se ||, and (b) V 1 x 10~ 3 (^/ 43 = 1.1). 



When a white Gaussian noise {SNR = 50dB) is added to both real and imaginary part of 
E and E l , the corresponding normalized vector is denoted by v n oise- Both ||v no j se || and ||v|| are 
plotted in Fig. [7J The missing cylinder is imaged and rj 1 is retrieved with a relative error less than 
0.1. The main impact of the addition of a white Gaussian noise to the data is a decrease of h l . 

As for the (i,j')-th DFSC configuration, the major difficulty is in the recovery of modified 
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cylinders close to one another. pfaw is symmetric in terms of q and I. Let us define the vector 
Q q such that Q q = Q l = mean/ (V^ q ' 1 ^, which enables us to depict the results in the same form as 
in the case of the (i)-th DFSC configuration. Figure E] shows both p(«>9 and Q q for a (43, 44)-th 
DFSC configuration with 9 = tt/2 when d = 2. Two identically modified cylinders are retrieved, 
separated from a distance d which is smaller than A/2. 



(a) (b) 




Figure 8: Crystal with two defects close to one another when d = 2: R = 0.15, rj = 1 and 
9 inc = vr/2. (a) x 1(T 3 (rf = 1.1) and (b) Q l . 

5 Conclusion 

The simple tomographic algorithm described herein appears as a first step in the monitoring of 
structured materials like PCs. Use of the specific Green's function together with the first-order 
Born approximation enables us to localize and characterize simple defects consisting in absence 
of cylinders or identical modification of the optical index of cylinders in finite-size PCs. Several 
results exhibit a possible retrieval of two defects beyond the Rayleigh criterion. The algorithm can 
be used as a first iteration in an iterative solution scheme. 

Various difficulties would arise if the low-frequency approximation was no more valid, which is 
the case for usual PCs at the location of band gap (their frequency band of use). The dependance 
of h! on the angle of incidence of the plane wave sollicitation as well as its dependance on the defect 
location is still a challenge. 
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